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Nanoparticles tethered with DNA strands are promising building blocks for bottom-up nanotech- 
nology, and a theoretical understanding is important for future development. Here we build on 
approaches developed in polymer physics to provide theoretical descriptions for the equilibrium 
clustering and dynamics, as well as the self-assembly kinetics of DNA-linked nanoparticles. Striking 
agreement is observed between the theory and molecular modeling of DNA tethered nanoparticles. 
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The specificity, directionality, and technological con- 
trol over base sequence make DNA an attractive linking 
unit for artificial constructs [TH5|. One such approach is 
to attach DNA strands of designed base sequence onto 
a nanoparticle (NP), thereby creating "molecules" that 
self-assemble into highly organized structures through 
complementary pairings of DNA [U |5]. Recent ex- 
periments have demonstrated that uniformly DNA- 
functionalized NP can self-assemble into disordered [BHS] 
or ordered crystal structures P^UHTB] . NPs with a discrete 
number of attachments are harder to prepare, but they 
provide possibilities for even more diverse structures |14l - 
I17j . In addition to examining the possible equilibrium 
structures, it is also vital to develop an understanding for 
the kinetics of the self-assembly process p. [T^ fT51 [TBI [T^ . 
since the pathway to desired structures may depend sen- 
sitively on sample preparation. 

In this Letter, we build on concepts from polymer 
physics to provide a theoretical framework that describes 
both the equilibrium properties and the self-assembly ki- 
netics of NPs functionalized with a small number of ss- 
DNA. We demonstrate the applicability of this theory 
to molecular simulations of a binary mixture where NPs 
are functionalized with two or three ssDNA, similar to 
experimentally realized systems [15j . 

The coarse-grained molecular model we use for the 
DNA-functionalized NP is a modest modification of a 
previous model [5D]. The inset of fig. [T] provides a graphi- 
cal representation. The sugar-phosphate backbone is rep- 
resented by connected beads with only excluded volume 
interactions; each bead carries an additional "sticky" site 
that represents a base (A, T, C, or G) and can bond only 
with another sticky site of complementary type (A-T or 
C-G pair). The size of the sticky site is small relative to 
the backbone beads, so that a base can bond to at most 
one other complementary base. A bending potential be- 
tween consecutive triplets of backbone beads is included 
to model the characteristic rigidity of the DNA strands 
and to maintain the angle between ssDNA attached to 
the same core (120° for 3- functional NP and 180° for 
2- functional NP). A more complete description of the 
model is provided in Ref. [20] . We consider ssDNA four 
bases in length with sequence A-C-G-T, chosen to en- 



able complementary pairing between ssDNA on different 
NPs. To mimic solvent effects on dynamics, we simulate 
this model using dissipative particle dynamics [21], an 
algorithm known to ensure correct hydrodynamic behav- 
ior. Length is measured in diameter of the beads, and 
temperature T is measured in units of the bonding en- 
ergy between complementary bases, taking Boltzmann's 
constant fc^ = 1 for simplicity. 

We consider a system with N^ = 945 2-functional par- 
ticles and TVs = 55 3-functional particles. The presence of 
3-functionalized particles facilitates branching of other- 
wise linear chains, thereby opening the possibility to form 
a percolating network and phase separation at low den- 
sity [22] . For part of our study where system size is less 
important, we consider a smaller system with N2 — 189 
and 7V3 = 11 so that we may more readily study the dy- 
namics of the most slowly relaxing systems. Nonetheless, 
equilibration at the lowest T studied in the small systenr 
still requires 5 x 10^ integration steps - roughly 3 months 
of computation on a single core with current resources. 
The number density is fixed at {N2+ N:i)/V = 0.0065 
(corresponding to mean separation of 5.4 bases between 
particles); we have confirmed with structural factor S{q) 
that there is no (gas-liquid) phase separation at this den- 
sity 

In this systenr, the extent of reaction p is characterized 
by the fraction of bonded strands. We define two strands 
as bonded when half or more of the base pairs on them are 
linked {i.e. the base-base potential energy is negative). 
The T dependence of equilibrium p is well described (see 
fig. [1] and Ref. [2D] ) by a two-state expression 
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where AU and AS" measure, respectively, the energy and 
entropy change associated with the formation of a double 
strand. 

To quantify the equilibrium assembly, we first calculate 
the number of clusters Nk of size k at each T (fig. [2la)). 
To describe our results, we consider the theoretical pre- 
dictions of Flory [23] and Stockmayer p4] (FS). The FS 
theory assumes that no intra-cluster loops are present 
and that every functional group is equally reactive. The 
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FIG. 1: Temperature dependence of the fraction of bonded 
arms at equilibrium, peq- Symbols are simulation data for 
system size N = 200 and A'^ = 1000. Solid line represents 
two-state expression with AU = 1.78 and AS' = 20.0. Dotted 
line indicates percolation threshold predicted by FS theory. 
Inset shows graphical representation of one 3-functional par- 
ticle bonding to one 2-functional particle. Large red spheres 
indicate NP, blue spheres indicate backbone beads of ssDNA, 
and small green spheres indicate location of "sticky" sites. 



theory predicts the most probable number of clusters Nni 
containing n 3-functional particles and / 2-functional par- 
ticles is 

Nniip) = mil -pr+'iP3pr-\p2pywni, 

^ {2n + l)\ (2) 



Wn 



n!n(n-f 2)!' 



where ps = "iNz/i^N^ + SiVa) and p2 = 1 — Pa are the 
portions of functional groups belonging to 3-functional or 
2-functional particles, respectively. A sum yields the an- 
alytic prediction Nk = X]n+i=fc -^"'- Since the only vari- 
able in eq. u^is the extent of reaction p (which we already 
know), we have a parameter free description of Nk- We 
find quantitative agreement between this prediction and 
our results (fig. [2]^a)). We can also evaluate the number 
of finite-size clusters Nc = X^fe-^fe; which we compare 
with the corresponding simulation results in fig. ^b). 
We again find quantitative agreement in the whole range 
of p - both below and above the percolation threshold 
Pc = 1/(1 + Ps) = 0.9256 - a demonstration that in this 
system an accurate analytic description for the equilib- 
rium clustering is possible. Qualitatively, these clusters 
are dominated by long chains of 2-functional NP, with 
occasional links provided by the 3-functional NP. 

The clustering process has a profound effect on the 
dynamics of the system, as revealed by the diffusion con- 
stant D, which drops several orders of magnitude in the 
small T-region where clusters form. Figure [3la) shows 
D evaluated from the asymptotic behavior of the mean- 
square displacement (Ar^) = 6Dt as a function of the 
inverse temperature 1/T. Close to the percolation tran- 
sition, D takes on an Arrhenius T dependence. 

To understand the behavior of D, consider that the 
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FIG. 2: (a) Equilibrium distribution of finite-size clusters A^^. 
Points are simulation data (system size A'' = 1000) and lines 
are the FS predictions. Dashed line represents the expected 
power-law decay Nk ~ k~^'^ at percolation [25]. (b) Num- 
ber of finite size clusters Nc as a function of p. Points are 
simulation data (system size A*' — 200) and black solid line 
is FS prediction. Inset enlarges the region above percolation 
threshold (indicated by dotted line). 



overall diffusion constant is an average over particles 
belonging to clusters of different sizes; therefore D — 
J2k^k{kNk)/N, where Dk is the diffusion constant for 
particles in a cluster of size k [26] . Typically, Dk drops 
inversely to the cluster size, so that D/, ^ Di/k, the 
so-called Stokesian limit of diffusion [?7l E5] . In this ap- 
proximation, D ^ Di J2k ^k/N = DiNc/N. If we take 
the monomer diffusion coefficient Di as a constant, we 
can readily test this prediction. Figure [3Fb) shows that 
D obeys this prediction well, with small deviations at 
the largest and smallest values. Numerically, we find 
Di = 0.145, roughly equal to the high T asymptotic 
value of D, where monomers dominate. As shown in 
fig. ^h), Nc ^ (1 — p) for p < Pc, and so we also expect 
D ^ (1 — p) to leading order (fig. Isla) inset). Since at 
low T, (1 — p) ^ expi—AU/ksT) (from eq. [I]), we can 
explain the Arrhenius behavior D ^ exp(— AL//A:bT) in 
the low T region. Indeed, AU obtained from p(T) (fi-g.fl]) 
matches that for D, demonstrating consistency. 

Having successfully described the clustering and dy- 
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FIG. 3: (a) Diffusion constant _D as a function of inverse tem- 
perature. The dotted line indicates the location of the perco- 
lation threshold, and red solid line represents Arrhenius be- 
havior with an activation energy AI7 — 1.78 - the same value 
reported for p{T) in fig.Il] See text for an explanation. Then 
inset shows D as a function of fraction of non-bonded arms, 
(b) Diffusion constant D versus number of finite clusters Nc- 
The line shows the predicted behavior D = Di Nc/N, where 
Di = 0.145. Calculation of (Ar^) is averaged over 8 to 15 
independent runs. System size A*' — 200. 



namics at equilibrium, we next explore to what degree 
the self-assembly kinetics can be analytically predicted. 
We focus on the evolution of the system after a tempera- 
ture jump (at t = 0) from a high-T unassociated state to 
a low-T self-assembled state, following the evolution of p 
from its initial value p(t = 0) « to its final equilibrium 
value p{t = oo) = Peq{T). With the assumption that ev- 
ery strand is equally reactive independent of the size of 
cluster to which it attaches, p{t) satisfies [5S] 
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where fcovcraii is the overall rate coefficient for bond break- 
ing. Equation [3] simply states that the rate of bond for- 
mation (first term) only depends on the probability of 
finding two unassociated strands, and that the rate of 
bond fragmentation (second term) only depends on the 
probability of finding a bonded strand. When this kinetic 



process is dominated by the intrinsic rate of reaction {i.e. 
dsDNA formation), fcovcraii = fcc is a constant, and eq. p^ 
can be integrated analytically to give the time evolution 
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where T — fcc(l + Pcq)/{1 — Pcq) determines the rate of 
reaction, and A = (1 ~ p{0)/pcq)/{l — p(O)pcq) incorpo- 
rates the initial condition p{0). It has been shown by van 
Dongen and Ernst [29 that, in this so-called "reaction- 
limited" or "chemical-limited" process, eq. [3] provides 
a time dependent distribution Nk{p(t)) (together with 
eq. ^ that satisfies the Smoluchowski rate equation with 
condensation and fragmentation terms below percola- 
tion I3(L 31J. 
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FIG. 4: (a) Rate of bond dp/dt as a function of fraction p of 
bonded strands, and (b) time evolution of p starting from a 
high-T state (Tinit = 0.2, p(0) = 0.008). Inset shows the two 
rate constants from best fit. Symbols are simulation results 
(system size N — 1000), averaged over 35 independent runs at 
short time and 15 runs at longer time. Dashed and solid lines 
represent, respectively, prediction of eq.|3]in the reaction-limit 
scenario (constant k) and with consideration of the pko rate 
coefficient. 

Using our simulation results, we test if our system can 
be described by the reaction-limited description given by 
eq.lSlwith fcovcraii = fcc- In fig-|4l we see that the reaction- 
limited expression describes the rate of bond dp/dt and 



time evolution p{t) well for most values of p but fails when 
p < 0.3. This indicates that, in the assembly process of 
DNA-functionalized NP, there is a rate-controlling factor 
for small p not captured by a simple reaction-limited de- 
scription. The trend in dp/dt shows that this additional 
rate-controlling factor slows down the overall reaction 
rate at small p. To incorporate the effect, we introduce 
an additional reaction rate /csmaii = ^oP that changes 
linearly with p. Since reaction rates sum inversely |27j . 
^overall ~ ^c^ + (^oP)~^- Figure. H shows that this modi- 
fied fcoveraii accurately describes our simulation results for 
the assembly kinetics over the entire self-assembly pro- 
cess. By numerically integrating eq. [3] using the appro- 
priate initial condition p{0) , we also provide an accurate 
estimate of the time evolution p{t) based on this expres- 
sion. Surprisingly, this prediction works well even for the 
case in which the final state percolates (T = 0.06). The 
inset in fig. |4] shows that both rate constants kc and fco 
follow an Arrhcnius behavior {i.e. Infc ~ 1/T). 

The presence of a second rate constant that is domi- 
nant at small p may be a consequence of the fact that 
bonding NP via dsDNA requires two steps: (i) the 
strands of opposing NP must be properly oriented, and 
(ii) once oriented, the linking of individual bases must 
occur. Since base-pair bonds are dynamic, there is a 
constant "flickering" of bonds. Thus, when p is signifi- 
cant (here, apparently p > 0.3), the flickering should be 
a consequence of the regular breaking and reforming of 
individual base-pair bonds. As a result, the strands will 
typically remain oriented to facilitate base-pair bonding, 
making the rate constant associated with orientation ir- 
relevant {i.e. fcovoraii dominated by k^). For very small 
p, the intermittency of the bonds is large enough that 
strands can "wander", and as a result the time needed 
to orient two strands becomes the rate limiting factor, 
and so fegmaii dominates. This is consistent with the 
fact that for simpler systems, where arm orientation is 
not a factor, there is only a single rate constant to con- 
sider [32] • The effect of this mechanism is the opposite of 
the diffusion-limited scenario, which is usually important 
only for long-time dynamics, when p is large [271 I30j. 

Our studies show that, by combining a thermodynamic 
description for Peq{T), the FS theory for Ni.{p), an ap- 
proximate expression for D, and rate equation with ap- 
propriate rate constants, we can provide a comprehensive 
theoretical description for the clustering, equilibrium dy- 
namics, and assembly kinetics that accurately predicts 
the behavior of this DNA-functionalized NP system. The 
complex process of self-assembly through dsDNA forma- 
tion can be simplified as a chemical-limited reaction plus 
an additional reaction rate kop associated with strand 
orientation. An obvious question is whether the results 
obtained for our simplified DNA simulation model will 
be transferable to real systems. The experimentally syn- 
thesized systems with very small numbers of attached 
DNA strands to a core NP studied by Alivisatos and 



co-workers fJS' offer the possibility to directly test our 
predictions. Additionally, our approach could be applica- 
ble to the experiments of Gang and co-workers ^ where 
there are many strands attached to a core NP, but only a 
small number of which are available for binding, resulting 
in more elongated clusters. They also consider micron- 
sized colloidal cores in that work, which can potentially 
significantly alter the assembly dynamics. 

The successful theoretical description of our numerical 
results shows that existing frameworks can be adapted 
describe novel self-assembled systems. The next step in 
this development is the satisfactory inclusion of closed 
loop structures, which is important for applications when 
the number of functionalizing strands is larger. 
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